/*
 * main.c
 * Main file for prime counting approximation code.  The goal is to obtain the most accurate approximation to pi(x) given limited time.
 * A formulation of Riemann's explicit formula provided by Tomas Oliveira e Silva is used, as well as a modification to it that
 * takes advantage of a statistical trick that involves average errors and a sieve of numbers around x.  Series expansions and complex analysis are used.
 *
 * Wikipedia lists prime counts for powers of 10 (August 2014):
 * x	pi(x)
 * 10	4
 * 10^2	25
 * 10^3	168
 * 10^4	1,229
 * 10^5	9,592
 * 10^6	78,498
 * 10^7	664,579
 * 10^8	5,761,455
 * 10^9	50,847,534
 * 10^10	455,052,511
 * 10^11	4,118,054,813
 * 10^12	37,607,912,018
 * 10^13	346,065,536,839
 * 10^14	3,204,941,750,802
 * 10^15	29,844,570,422,669
 * 10^16	279,238,341,033,925
 * 10^17	2,623,557,157,654,233
 * 10^18	24,739,954,287,740,860
 * 10^19	234,057,667,276,344,607
 * 10^20	2,220,819,602,560,918,840
 * 10^21	21,127,269,486,018,731,928
 * 10^22	201,467,286,689,315,906,290
 * 10^23	1,925,320,391,606,803,968,923
 * 10^24	18,435,599,767,349,200,867,866
 * 10^25	176,846,309,399,143,769,411,680
 *
 * Two approximations are considered here.  The first is a formulation of Riemann's explicit formula provided by Tomas Oliveira e Silva, and the second
 * is a modification to Riemann's explicit formula which sieves an interval around x and tabulates an average to arrive at a formula that can be
 * expected to converge more quickly than Riemann's explicit formula alone.
 * Read the documentation of the functions for more information.
 * Also read Kevin Stueve's Winter 2010 Number Theory paper "Monte-Carlo Approximation of the Prime Counting Function."
 *      Author: Kevin Marshall Stueve; August 2014
 * Copyright (C) 2014 Kevin Marshall Stueve
 *
 * This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <limits.h>
#include "real_li.h"
#include "complex_li.h"

double prime_powers_correction(long long x,long long lower_radius,long long upper_radius);
int* generate_primes_up_to(int x, int* number_primes);
long double li_int(long double x);
char* create_initial_sieving_array();
int sum_deltas(long long x,int lower_radius,int upper_radius);
void sieve(char sieving_array[],int sieving_array_size,int starting_x);
void sieve_for_prime_power_correction(char sieving_array[],int sieving_array_size);
void tear_down();
void riemann_init(void);
void riemann_test();
void pi_approx_init();
void pi_approx_test();
char* initial_sieving_array;
int initial_sieving_array_length;
int* coprime_integers;
int coprime_integers_length;
int wheel_size;
int first_sieving_prime;
int* sieving_primes;
int number_sieving_primes;

/*
 * Main entry point into the code.
 * real_li_test
 *     li_approx
 *     li
 * riemann_init
 *     import_zeros
 * riemann_pi
 *     li_rho
 *     simple_pi
 * pi_approx_init
 *     generate_primes_up_to
 *     create_initial_sieving_array
 * pi_approx
 *     li_int
 *     x_li_rho
 *     li_rho_plus_one
 *     prime_powers_correction
 *         sieve_for_prime_power_correction
 *     sum_deltas
 *         sieve
 *
 *     Author: Kevin Marshall Stueve; August 2014
 */
int main(void) {
	real_li_test();
	printf("\n");

	riemann_init();
	riemann_test();
    printf("\n");

	pi_approx_init();
	pi_approx_test();
	printf("\n");

    printf("simple_pi(10^6+1)-simple_pi(10^6):%d\n",simple_pi(1000001)-simple_pi(1000000));
    printf("simple_pi(10^9+1)-simple_pi(10^9):%d\n",simple_pi(1000000001)-simple_pi(1000000000));
    printf("\n");

    tear_down();

	printf("end");
}

/*
 * Release memory back to the OS
 *     Author: Kevin Marshall Stueve; August 2014
 */
void tear_down() {
	free(sieving_primes);
	free(coprime_integers);
}

/*
 * A formulation of Riemann's explicit formula for computerized calculation where the influence of prime powers is absorbed into recursive pi(x) calls.
 * From Tomas Oliveira e Silva:
 * "pi(x)=li(x)
 *      -sum_{\rho} li(x^{\rho})
 *      -log(2)
 *      -sum_{k=2}^{\infty} pi(x^{1/k})/k
 *      +\int_x^\infty dt/(t*(t^2-1)*log(t))"
 *
 * This is equivalent to Riemann's explicit formula.
 * li(x) is the logarithmic integral (starting at 0)
 * rho is taken to be 1/2 +/- it, where t runs over the nontrivial zeros of zeta, the Riemann zeta function
 *     (the formula would still be valid if a rho was off the critical line real=1/2)
 * log is taken here to be the natural logarithm, with base e
 * The sum in the 4th line (with the recursive pi(x) calls) terminates after a finite number of terms.
 * The final term is included for correctness, but can be removed in the final calculation because of its small size.
 */
long double riemann_pi(long long int x,int num_zeros) {
	long double result = li(x) - log(2);

	for (int i=0;i<num_zeros;i++){
		result -= li_rho(x,zeta_zeros[i]);
	}
	double k=2;
	long double term;
	do {
      term=simple_pi(pow(x,1.0/k))/k;
      result-=term;
      k++;
	} while (term>0);
	return result;
}
/*
 * Simple pi(x), for the recursive pi(x) calls in the above formula, riemann_pi (unoptimized)
 *     Author: Kevin Marshall Stueve; August 2014
 */
long long int simple_pi(long long int x) {
    //int sieving_array[x>>1];

    char * sieving_array = (char *) malloc(x>>1);

    for (int i=1; i < x; i+=2) {
    	sieving_array[i>>1] = 1;
    }
    int sieving_prime = 3;
    int sieving_prime_limit = floor(sqrt(x));
    while (sieving_prime<=sieving_prime_limit) {
    	for (int i=sieving_prime*sieving_prime; i<x; i += 2*sieving_prime) {
    		sieving_array[i>>1] = 0;
    	}
        do {
        	sieving_prime+=2;
        } while (sieving_array[sieving_prime>>1]==0);
    }
    int count = 0;
    for (int i = 1; i < x; i+=2) {
    	count += sieving_array[i>>1];
    }
    free(sieving_array);
    return count;
}


/*
 * Before any calculation, we must import zeros of the Riemann zeta function into the program.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void riemann_init(void) {
	NUM_ZEROS=import_zeros();
}

/*
 * An array of tests of the riemann_pi function.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void riemann_test() {
	printf("riemann_pi test\n");
	printf("num_zeros:%d\n",NUM_ZEROS);
	printf("zero 0: %f\n",zeta_zeros[0]);
	printf("zero 2001051: %f\n",zeta_zeros[2001051]);
    printf("\n");

	printf("simple_pi(100): %llu\n",simple_pi(100));
	printf("riemann_pi(100): %Lf\n", riemann_pi(100,NUM_ZEROS));
	printf("\n");

	printf("simple_pi(1000): %llu\n",simple_pi(1000));
	printf("riemann_pi(1000): %Lf\n", riemann_pi(1000,NUM_ZEROS));
	printf("\n");

	printf("simple_pi(10000): %llu\n",simple_pi(10000));
	printf("riemann_pi(10000): %Lf\n", riemann_pi(10000,NUM_ZEROS));
	printf("\n");

	printf("simple_pi(100000): %llu\n",simple_pi(100000));
	printf("riemann_pi(100000): %Lf\n", riemann_pi(100000,NUM_ZEROS));
	printf("\n");

	printf("simple_pi(1000000): %llu\n",simple_pi(1000000));
	printf("riemann_pi(1000000): %Lf\n", riemann_pi(1000000,NUM_ZEROS));
	printf("\n");

	//printf("simple_pi(10000000): %f\n",simple_pi(10000000));
	//printf("riemann_pi(10000000): %f\n", riemann_pi(10000000));

	//printf("simple_pi(100000000): %f\n",simple_pi(100000000));
	//printf("riemann_pi(100000000): %f\n", riemann_pi(100000000));

	//printf("simple_pi(1000000000): %f\n",simple_pi(1000000000));
	//printf("riemann_pi(1000000000): %f\n", riemann_pi(1000000000));

}


/*
 * An approximation to pi(x) to be studied discovered by Kevin Marshall Stueve.  Compare to riemeann_pi.
 *     Author: Kevin Marshall Stueve; August 2014
 */long double pi_approx(long long x,long long num_zeros,long long lower_radius,long long upper_radius) {
	 if (x<0 || num_zeros<0 || lower_radius<0 ||upper_radius <0) {
	 	printf("Error in average_deltas, inputs must be positive.\n");
	 	return 0;
	 }
	 //The starting point is the average of the logarithmic integral.
	 long double result = (li_int(x+upper_radius+1)
			             -li_int(x-lower_radius)) / (upper_radius+lower_radius+1);

	//We must include the -log(2) term from Riemann's explicit formula
	result-=log(2);

	//Int li(x^rho)= x*li(x^rho)+ li(x^(rho+1))
	long double sum=0;
	for (int i=0;i<num_zeros;i++){
		sum+=(x_li_rho(x+upper_radius+1,zeta_zeros[i])
			 -x_li_rho(x-lower_radius,  zeta_zeros[i])) / (upper_radius+lower_radius+1);
    }
	result-=sum;
	sum=0;
	for (int i=0;i<num_zeros;i++){
		sum+=(li_rho_plus_one(x+upper_radius+1,zeta_zeros[i])
			 -li_rho_plus_one(x-lower_radius,  zeta_zeros[i])) / (upper_radius+lower_radius+1);
	}
    result+=sum;

	//result-=Avg sum_{k=2}^{\infty} pi(x^{1/k}/k)
	result-=prime_powers_correction(x,lower_radius,upper_radius) / (upper_radius+lower_radius+1);

	//result+=Avg pi(x)-pi(x_i)
	result+=sum_deltas(x,lower_radius,upper_radius) / (upper_radius+lower_radius+1);

	return result;
}
 /*Steps which must be performed prior to calling pi_approx.
  *    Author: Kevin Marshall Stueve; August 2014
  */
void pi_approx_init() {
	sieving_primes = generate_primes_up_to(sqrt(1e9+1e9),(int*)&number_sieving_primes);//1e9+1e9
	initial_sieving_array = create_initial_sieving_array();
}

/* A comparison of riemann_pi to pi_approx for x=1000000 and n=1000.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_1(){
	printf("pi_approx_test_1, 1000 zeros\n");
	printf("riemann_pi(1000000,10000):             %Lf\n",riemann_pi(1000000,1000)-78498);
	printf("pi_approx(1000000,1000,1,1):           %Lf\n", pi_approx(1000000,1000,1,1)-78498);
	printf("pi_approx(1000000,1000,10,10):         %Lf\n", pi_approx(1000000,1000,10,10)-78498);
	printf("pi_approx(1000000,1000,100,100):       %Lf\n", pi_approx(1000000,1000,100,100)-78498);
	printf("pi_approx(1000000,1000,1000,1000):     %Lf\n", pi_approx(1000000,1000,1000,1000)-78498);
	printf("pi_approx(1000000,1000,10000,10000):   %Lf\n", pi_approx(1000000,1000,10000,10000)-78498);
	printf("pi_approx(1000000,1000,100000,100000): %Lf\n", pi_approx(1000000,1000,100000,100000)-78498);
}

/* A comparison of riemann_pi to pi_approx for x=1000000 and n=10000.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_2(){
	printf("pi_approx_test_2, 10000 zeros\n");
	printf("riemann_pi(1000000,10000):              %Lf\n",riemann_pi(1000000,10000)-78498);
	printf("pi_approx(1000000,10000,1,1):           %Lf\n", pi_approx(1000000,10000,1,1)-78498);
	printf("pi_approx(1000000,10000,10,10):         %Lf\n", pi_approx(1000000,10000,10,10)-78498);
	printf("pi_approx(1000000,10000,100,100):       %Lf\n", pi_approx(1000000,10000,100,100)-78498);
	printf("pi_approx(1000000,10000,1000,1000):     %Lf\n", pi_approx(1000000,10000,1000,1000)-78498);
	printf("pi_approx(1000000,10000,10000,10000):   %Lf\n", pi_approx(1000000,10000,10000,10000)-78498);
	printf("pi_approx(1000000,10000,100000,100000): %Lf\n", pi_approx(1000000,10000,100000,100000)-78498);
}
/* A comparison of riemann_pi to pi_approx for x=1000000 and n=100000.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_3(){
	printf("pi_approx_test_3, 100000 zeros\n");
	printf("riemann_pi(1000000,100000):              %Lf\n",riemann_pi(1000000,100000)-78498);
	printf("pi_approx(1000000,100000,1,1):           %Lf\n", pi_approx(1000000,100000,1,1)-78498);
	printf("pi_approx(1000000,100000,10,10):         %Lf\n", pi_approx(1000000,100000,10,10)-78498);
	printf("pi_approx(1000000,100000,100,100):       %Lf\n", pi_approx(1000000,100000,100,100)-78498);
	printf("pi_approx(1000000,100000,1000,1000):     %Lf\n", pi_approx(1000000,100000,1000,1000)-78498);
	printf("pi_approx(1000000,100000,10000,10000):   %Lf\n", pi_approx(1000000,100000,10000,10000)-78498);
	printf("pi_approx(1000000,100000,100000,100000): %Lf\n", pi_approx(1000000,100000,100000,100000)-78498);
}

/* A comparison of riemann_pi to pi_approx for x=1000000 and n=1000000.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_4(){
	printf("pi_approx_test_4, 1000000 zeros\n");
	printf("riemann_pi(1000000,1000000):              %Lf\n",riemann_pi(1000000,1000000)-78498);
	printf("pi_approx(1000000,1000000,1,1):           %Lf\n", pi_approx(1000000,1000000,1,1)-78498);
	printf("pi_approx(1000000,1000000,10,10):         %Lf\n", pi_approx(1000000,1000000,10,10)-78498);
	printf("pi_approx(1000000,1000000,100,100):       %Lf\n", pi_approx(1000000,1000000,100,100)-78498);
	printf("pi_approx(1000000,1000000,1000,1000):     %Lf\n", pi_approx(1000000,1000000,1000,1000)-78498);
	printf("pi_approx(1000000,1000000,10000,10000):   %Lf\n", pi_approx(1000000,1000000,10000,10000)-78498);
	printf("pi_approx(1000000,1000000,100000,100000): %Lf\n", pi_approx(1000000,1000000,100000,100000)-78498);
}

/* A comparison of riemann_pi to pi_approx for x=1000000 and n=1000000.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_5(){
	printf("pi_approx_test_5,1000000 zeros\n");
	printf("riemann_pi(1000000,1000000):    %Lf\n",riemann_pi(1000000,1000000)-78498);
	printf("pi_approx(1000000,1000000,1,1): %Lf\n", pi_approx(1000000,1000000,1,1)-78498);
	printf("pi_approx(1000000,1000000,0,0): %Lf\n", pi_approx(1000000,1000000,0,0)-78498);
	printf("pi_approx(1000000,1000000,1,0): %Lf\n", pi_approx(1000000,1000000,1,0)-78498);
	printf("pi_approx(1000000,1000000,0,1): %Lf\n", pi_approx(1000000,1000000,0,1)-78498);
}

/* A comparison of riemann_pi to pi_approx for x=1000000000 and n=2001052.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_6() {
	printf("pi_approx_test_6, 2001052 zeros\n");
	printf("riemann_pi(1000000000,2001052):    %Lf\n",riemann_pi(1000000000,2001052)-50847534);
	printf("pi_approx(1000000000,2001052,1,1): %Lf\n", pi_approx(1000000000,2001052,1,1)-50847534);
	printf("pi_approx(1000000000,2001052,0,0): %Lf\n", pi_approx(1000000000,2001052,0,0)-50847534);
	printf("pi_approx(1000000000,2001052,1,0): %Lf\n", pi_approx(1000000000,2001052,1,0)-50847534);
	printf("pi_approx(1000000000,2001052,0,1): %Lf\n", pi_approx(1000000000,2001052,0,1)-50847534);
}

/* A comparison of riemann_pi to pi_approx for x=1000000000 and n=2001052.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_7() {
	printf("pi_approx_test_7\n");
	printf("riemann_pi(1000000000,2001052):                    %Lf\n",riemann_pi(1000000000,2001052)-50847534);
	printf("pi_approx(1000000000,2001052,10,10):               %Lf\n", pi_approx(1000000000,2001052,10,10)-50847534);
	printf("pi_approx(1000000000,2001052,100,100):             %Lf\n", pi_approx(1000000000,2001052,100,100)-50847534);
	printf("pi_approx(1000000000,2001052,1000,1000):           %Lf\n", pi_approx(1000000000,2001052,1000,1000)-50847534);
	printf("pi_approx(1000000000,2001052,10000,10000):         %Lf\n", pi_approx(1000000000,2001052,10000,10000)-50847534);
	printf("pi_approx(1000000000,2001052,100000,100000):       %Lf\n", pi_approx(1000000000,2001052,100000,100000)-50847534);
	printf("pi_approx(1000000000,2001052,1000000,1000000):     %Lf\n", pi_approx(1000000000,2001052,1000000,1000000)-50847534);
	printf("pi_approx(1000000000,2001052,10000000,10000000):   %Lf\n", pi_approx(1000000000,2001052,10000000,10000000)-50847534);
	printf("pi_approx(1000000000,2001052,100000000,100000000): %Lf\n", pi_approx(1000000000,2001052,100000000,100000000)-50847534);
}

/* A comparison of riemann_pi to pi_approx for x=1000000000000 and n=2001052.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test_8() {
	printf("pi_approx_test_8\n");
	printf("riemann_pi(1000000000000,2001052):    %Lf\n",riemann_pi(1000000000000,2001052)-37607912018);
	printf("pi_approx(1000000000000,2001052,1,1): %Lf\n", pi_approx(1000000000000,2001052,1,1)-37607912018);
	printf("pi_approx(1000000000000,2001052,0,0): %Lf\n", pi_approx(1000000000000,2001052,0,0)-37607912018);
	printf("pi_approx(1000000000000,2001052,1,0): %Lf\n", pi_approx(1000000000000,2001052,1,0)-37607912018);
	printf("pi_approx(1000000000000,2001052,0,1): %Lf\n", pi_approx(1000000000000,2001052,0,1)-37607912018);
}

/*An assortment of comparisons between pi_approx and riemann_pi
 *    Author: Kevin Marshall Stueve; August 2014
 */
void pi_approx_test() {
	printf("Errors:\n");
	pi_approx_test_1();
	printf("\n");
   pi_approx_test_2();
	printf("\n");
	pi_approx_test_3();
    printf("\n");
    pi_approx_test_4();
    printf("\n");
    pi_approx_test_5();
    printf("\n");
    pi_approx_test_6();
    //printf("\n");
    //pi_approx_test_7();
    //printf("\n");
    //pi_approx_test_8();

    //printf("riemann_pi(10^8):%f\n",riemann_pi(pow(10,8),1000));
    //printf("pi_approx(10^8) :%f\n", pi_approx(pow(10,8),2001052,pow(10,6),pow(10,6)));
	//printf("riemann_pi(10^9):%f\n",riemann_pi(pow(10,9),1000));
    //printf("pi_approx(10^9) :%f\n", pi_approx(pow(10,9),2001052,pow(10,6),pow(10,6)));

    //printf("riemann_pi(%u):%f\n",10000000,riemann_pi(10000000));
	//printf("pi_approx(%u) :%f\n",10000000, pi_approx(10000000,2001052,1000,1000));//10,10  100 100

	//printf("pi_approx:%f\n",pi_approx(10000,2001052,10,10));//100,100
    //printf("pi_approx:%f\n",pi_approx(1000000,2001052,10,10));//100,100
    //printf("pi_approx:%f\n",pi_approx(1000000,2001052,100,100));
	//printf("pi_approx:%f\n",pi_approx(1000000,2001052,10000,10000));
	//printf("pi_approx:%f\n",pi_approx(1000000,2001052,100000,100000));
}

/*
 * Perform an initial sieving by primes up to 23 in a block of size 510510 that can be used repeatedly
 *     Author: Kevin Marshall Stueve; August 2014
 */
 char* create_initial_sieving_array() {
 	int primes[]={2,3,5,7,11,13,17,19,23};
 	wheel_size = 2*3*5*7*11*13*17*19*23; //510,510
    first_sieving_prime = 29;

 	coprime_integers_length=(2-1)*(3-1)*(5-1)*(7-1)*(11-1)*(13-1)*(17-1)*(19-1)*(23-1);
    coprime_integers = (int *) malloc(coprime_integers_length * sizeof(int));

 	char * result = (char *) malloc(wheel_size);
 	initial_sieving_array_length = wheel_size;

 	//initialize array to 1's
 	for (int i=0;i<wheel_size;i++){
 		result[i]=1;
 	}

 	result[0]=0;//0 is not counted as prime
 	//Note that we must remember to set result[1]=0 when using the result of this function near the origin

 	//sieve by primes
 	for (int i=0;i<sizeof(primes)/sizeof(primes[0]);i++) {
 		for (int j=primes[i];j<wheel_size;j+=primes[i]) {//inefficiency
 			result[j]=0;
 		}
 	}
 	int coprime_integers_index=0;
 	for (int i=0; i<wheel_size;i++) {
 		if (result[i]==1){
 			coprime_integers[coprime_integers_index]=i;
 			coprime_integers_index++;
 		}
 	}
 	return result;
 }

 /*
  * Calculates the contribution of prime powers to Riemann's explicit formula in pi_approx.
  * That is -Avg sum_{k=2}^{\infty} pi(x^{1/k})/k
  *     Author: Kevin Marshall Stueve
  */
 double prime_powers_correction(long long x,long long lower_radius,long long upper_radius) {
 	 int v = sqrt(x+upper_radius)+1;
 	 int array_size = v;
 	 char * sieving_array = (char *) malloc(array_size);
 	 sieve_for_prime_power_correction(sieving_array,array_size);

 	 int counts [array_size];
 	 int count =0;

 	 for (int i=0;i<=v;i++) {
 		 count+=sieving_array[i];
 		 counts[i] = count;
 	 }

 	 long double sum = 0;
 	 for (int k=2;x+upper_radius>=pow(2,k);k++) {
          for(int i = -lower_radius;i<=upper_radius;i++ ) {
         	 v=pow(x+i,1.0/k);
         	 sum+=counts[v]/(k+0.0);
          }
 	 }
 	 free(sieving_array);
 	 return sum;
  }

 /*
   * To compute pi(x), we will call generate_primes_up_to(sqrt(x+upper_limit),&number_primes) to obtain primes up to sqrt(x+upper_radius).
   * Within generate_primes_up_to_x, primes up to fourth_root(x+upper_radius) will be used to find primes up to sqrt(x+_upper_radius)
   *     Author: Kevin Marshall Stueve; August 2014
   */
int* generate_primes_up_to(int x, int* number_primes) {
	int num_primes_up_to_x = simple_pi(x);

	//result is a list of primes up to x
	int *result = malloc(num_primes_up_to_x * sizeof(int));

	int sieving_limit = x;
	int sieving_prime_limit = floor(sqrt(x));

	//initialize sieving_array
	char *sieving_array = malloc(sieving_limit + 1);
	for (int i = 0; i <= sieving_limit; i++) {
		sieving_array[i] = 1;
	}

	//sieve out 0 and 1
	sieving_array[0] = 0;
	sieving_array[1] = 0;

	//sieve out composites
	for (int sieving_prime = 2; sieving_prime <= sieving_prime_limit; sieving_prime += 1) { 	//inefficient
		for (int i = sieving_prime * sieving_prime; i <= sieving_limit; i +=sieving_prime) {
			sieving_array[i] = 0;
		}
	}

	//compile primes into array
	int primes_index = 0;
	for (int i = 0; i <= sieving_limit; i++) {
		if (sieving_array[i] == 1) {
			result[primes_index] = i;
			primes_index += 1;
		}
	}

	//update number of primes in array
	*number_primes = primes_index;

	free(sieving_array);
	return result;
}

 /*
  * Integral of li(x).
  *    Author: Kevin Marshall Stueve; August 2014
  */
 long double li_int(long double x) {
 	return x*li(x)-li(x*x);
 }

 /*
  * The maximum of two numbers.
  *     Author: Kevin Marshall Stueve; August 2014
  */
long double max(long double x, long double y){
	if (x>y) {
		return x;
	}
	return y;
}

/*
 * Given a sieving array of some size, where the start of the array represents some number, sieve composites from the array.
 * Used by sum_deltas.
 *     Author: Kevin Marshall Stueve; August 2014
 */
void sieve(char sieving_array[],int sieving_array_size,int starting_x) {
     for (int i =0; i<sieving_array_size;i++){
     	sieving_array[i]=initial_sieving_array[(i+starting_x)%initial_sieving_array_length];
     }

     int sieving_primes_index = 0;
     while (sieving_primes[sieving_primes_index]!=first_sieving_prime) {
     	sieving_primes_index++;
     }

     for (;sieving_primes[sieving_primes_index]*sieving_primes[sieving_primes_index]<starting_x+sieving_array_size;sieving_primes_index++) {
     	int i=starting_x-starting_x%sieving_primes[sieving_primes_index];
     	if (i<starting_x) {
     		i+=sieving_primes[sieving_primes_index];
     	}
     	i = max(i,sieving_primes[sieving_primes_index]*sieving_primes[sieving_primes_index]);
     	while (i<starting_x+sieving_array_size) {
 			sieving_array[i-starting_x]=0;
 			i+=sieving_primes[sieving_primes_index];
     	}
     }
 }

 /*
  * Given a sieving array of some size, starting at the origin, sieve composites from the array.
  * Used by prime_powers_correction.
  *     Author: Kevin Marshall Stueve; August 2014
  */
 void sieve_for_prime_power_correction(char sieving_array[],int sieving_array_size) {
      for (int i =0; i<sieving_array_size;i++){
      	sieving_array[i]=initial_sieving_array[i%initial_sieving_array_length];
      }
      sieving_array[1]=0;
      if (sieving_array_size>2)
          sieving_array[2]=1;
      if (sieving_array_size>3)
          sieving_array[3]=1;
      if (sieving_array_size>5)
          sieving_array[5]=1;
      if (sieving_array_size>7)
          sieving_array[7]=1;
      if (sieving_array_size>=11)
          sieving_array[11]=1;
      if (sieving_array_size>13)
          sieving_array[13]=1;
      if (sieving_array_size>17)
          sieving_array[17]=1;
      if (sieving_array_size>19)
          sieving_array[19]=1;
      if (sieving_array_size>23)
          sieving_array[23]=1;

      int sieving_primes_index = 0;
      while (sieving_primes[sieving_primes_index] != first_sieving_prime) {
          sieving_primes_index++;
      }

      for (;sieving_primes[sieving_primes_index]*sieving_primes[sieving_primes_index]<sieving_array_size;sieving_primes_index++) {
      	int i=sieving_primes[sieving_primes_index]*sieving_primes[sieving_primes_index];
      	while (i<sieving_array_size) {
  			sieving_array[i]=0;
  			i+=sieving_primes[sieving_primes_index];
      	}
      }
 }

 /*
  * Sum(pi(x)-pi(x_i))
  * x_i ranges from x-lower_radius to x+upper_radius, inclusive on both ends.
  *     Author: Kevin Marshall Stueve; August 2014
  */
 int sum_deltas(long long x,int lower_radius,int upper_radius) {
 	int sieving_array_size = lower_radius+1+upper_radius;//+2
 	char * sieving_array = (char *) malloc(sieving_array_size);

 	sieve(sieving_array,sieving_array_size,x-lower_radius);

    int total_sum=0;
    int count =0;
 	for (int i=lower_radius;i>0;i--) {
 	    count+=sieving_array[i];
 	    total_sum+=count;
 	}
 	int result = total_sum;

 	total_sum=0;
 	count = 0;
 	for (int i = lower_radius+1;i<=lower_radius+upper_radius;i++){
 		count+=sieving_array[i];
 		total_sum-=count;
 	}
 	result+=total_sum;
 	free(sieving_array);
 	return result;
 }
